##########################################################################################

library(data.table)
library(optparse)
library(parallel)
library(dplyr)
library(RColorBrewer)
library(ggplot2)
library(ggrepel)

##########################################################################################

option_list <- list(
    make_option(c("--sample_list_file"), type = "character"),
    make_option(c("--rsem_file"), type = "character"),
    make_option(c("--gtf_file"), type = "character"),
    make_option(c("--out_path"), type = "character")
)

if(1!=1){
    
    sample_list_file <- "~/20220915_gastric_multiple/rna_batch1/analysis/config/tumor_normal.list"
    diff_file <- "~/20220915_gastric_multiple/rna_batch1/analysis/images/DiffGene/DiffGene.tsv"
    out_path <- "~/20220915_gastric_multiple/rna_batch1/analysis/images/Pathway"
    pathway_path <- "~/ref/Pathway/"
    gtf_file <- "~/ref/GTF/gencode.v19.ensg_genename.txt"

}

parseobj <- OptionParser(option_list=option_list, usage = "usage: Rscript %prog [options]")
opt <- parse_args(parseobj)
print(opt)

sample_list_file <- opt$sample_list_file
rsem_file <- opt$rsem_file
out_path <- opt$out_path
gtf_file <- opt$gtf_file

dir.create(out_path , recursive = T)

##########################################################################################

info <- data.frame(fread(sample_list_file))
dat_gtf <- data.frame(fread(gtf_file , header = F))
colnames(dat_gtf) <- c("gene_id" , "Hugo_Symbol")

result_wilcox <- data.frame(fread(diff_file , header = T))

##########################################################################################

class_type <- c( "Normal" , "IM" , "IGC" , "DGC")

##########################################################################################

diffplot <- function(use_dat = use_dat , image_name = image_name){

    #对原数据进行处理
    use_dat$padj <- -log10(use_dat$p)
    use_dat$log_FC <- log2(1/use_dat$foldchange)
    use_dat$gene <- use_dat$Hugo_Symbol
    use_dat$log_FC <- ifelse( use_dat$log_FC > 5 , 5 , use_dat$log_FC )
    use_dat$log_FC <- ifelse( use_dat$log_FC < -5 , -5 , use_dat$log_FC )

    data <- use_dat[,c('gene','log_FC','padj')]
    colnames(data) <- c('gene','log_FC','-log10_P_Value')

    #设置阈值
    logFC_cutoff <- log2(1.5)
    log10_P_Value_cutoff <- -log10(0.05)

    plot1 <- ggplot(data = data,aes(x = log_FC,y = `-log10_P_Value`))+
      geom_point(data = subset(data,abs(log_FC)<logFC_cutoff),
                 aes(size = abs(log_FC)),col = 'gray',alpha = 0.4)+
      geom_point(data = subset(data,abs(`-log10_P_Value`)<log10_P_Value_cutoff & abs(log_FC)>logFC_cutoff),
                 aes(size = abs(log_FC)),col = 'gray',alpha = 0.4)+
      geom_point(data = subset(data,abs(`-log10_P_Value`)>log10_P_Value_cutoff & log_FC>logFC_cutoff),
                 aes(size = abs(log_FC)),col = 'red',alpha = 0.4)+
      geom_point(data = subset(data,abs(`-log10_P_Value`)>log10_P_Value_cutoff & log_FC< -logFC_cutoff),
                 aes(size = abs(log_FC)),col = 'darkgreen',alpha = 0.4)+
      theme_bw()+
      theme(legend.title = element_blank(),
            panel.grid.major = element_blank(),
            panel.grid.minor = element_blank(),
            legend.position = 'none',
            axis.line = element_line(colour = "black"))+
      labs(x='log2(fold change)',y='-log10(adjusted p-value)')+
      geom_vline(xintercept = c(-logFC_cutoff,logFC_cutoff),lty = 3,col = 'black',lwd = 0.4)+
      geom_hline(yintercept = log10_P_Value_cutoff,lty = 3,col = 'black',lwd = 0.4) +
      geom_text_repel(data = subset(data,abs(`-log10_P_Value`)>log10_P_Value_cutoff & abs(log_FC)>logFC_cutoff),
                            aes(label = gene),size = 5,col = 'black')
    ggsave( image_name , plot1 )
}

##########################################################################################
## 分高表达和低表达，分布做通路富集
use_dat <- subset( result_wilcox , class1=="Normal" & class2=="IM")
image_name <- paste0( out_path , "/IM_Normal.DiffGene.pdf" )
diffplot(use_dat = use_dat , image_name = image_name)

use_dat <- subset( result_wilcox , class1=="Normal" & class2=="IGC")
image_name <- paste0( out_path , "/IGC_Normal.DiffGene.pdf" )
diffplot(use_dat = use_dat , image_name = image_name)

use_dat <- subset( result_wilcox , class1=="Normal" & class2=="DGC")
image_name <- paste0( out_path , "/DGC_Normal.DiffGene.pdf" )
diffplot(use_dat = use_dat , image_name = image_name)

